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1. Introduction 

To investigate the statistical properties of a many-body interacting system, often 
sampling following the predetermined distribntion is employed. Becanse of strong 
correlation between degrees of freedom, it is usnally intractable to compnte the 
expectation, covariance, etc. following the distribntion by direct manipnlation. Instead, 
we often use a sampling method. One of the well-known techniques is the Markov- 
chain Monte Carlo (MCMC) method, which imitates the stochastic dynamics of 
the discrete-time master equation. This is a promising method used to evaluate 
the desired quantities, but it sometimes takes a long time to convert the initially 
prepared system into a suitable state for computation because the method imitates 
the relaxation of the stochastic dynamics to the predetermined distribution as a steady 
state. Innovative studies to enhance its sampling efficacy and accelerate the relaxation 
to the predetermined distribution have been extensively reported niEiEiiaEiiniiTii. 

Basically, we construct a transition matrix to relax the initially prepared system 
to the steady state, which corresponds to the predetermined distribution in the 
MCMC method. A series of instantaneous realizations is generated by simulating the 
multiplication of the transition matrix. After sufficient iteration steps, the obtained 
series of realizations follows the predetermined distribution. This is assured by the 
balance condition (BC) on the transition matrix. In general, it is difficult to solve the 
BC to compute the transition matrix. The standard approach for hnding the solution is 
by use of the detailed balance condition (DBC). Because the DBC always satishes the 
BC, we may instead solve the DBC for construction of the transition matrix. However, 
the solution under the DBC is not necessarily an optimal choice of the transition matrix. 
Recent studies reveal the superior performance, which is conhrmed by means of the faster 
convergence to the predetermined distribution, by using nontrivial solution neglecting 
the DBC [nUinilll]. As detailed also in the present paper, the present authors provided 
its mathematical validation through an analysis on the eigenvalues of the transition 
matrix in a previous study [12] . 

The method for sampling is not limited to the MCMC method. The simplest 
technique is to implement Langevin dynamics, which is available for the stochastic 
dynamics of continuous variables in continuous-time evolution. In this alternative, 
does the analogous procedure for achieving faster convergence to the predetermined 
distribution exist? The time evolution of the probabilistic distribution for the degrees 
of freedom driven by Langevin dynamics is written by using the Fokker-Planck equation. 
This is the continuous equation for the instantaneous distribution function. In the steady 
state, the divergence of the flow of the distribution must vanish. This is analogous to 
the BC in the context of the master equation. The condition that the steady flow itself 
vanishes corresponds to the DBC. We often demand that the force, which drives the 
system to simulate the natural dynamics, be the gradient of the potential energy. This 
natural force holds the analogous condition to the BC and can be conhrmed to satisfy 
the DBC in a different formulation. However, it is not necessarily the case if we desire 


Mathematical understanding of violation of detailed balance condition 


3 


to reach the predetermined distribution as fast as possible. We have room to design a 
force to speed up the convergence to the desired distribution. In the present study, we 
propose several forces to violate the condition analogous to the DBG while satisfying 
the analogous condition to the BC in the steady state. 

The rest of the paper consists of the following sections. We review several ways 
to construct the transition matrix neglecting the DBG in the MGMG method in the 
next section. After that, we introduce the Langevin equation and the corresponding 
Fokker-Planck equation in the third section. In this section, we hnd several nontrivial 
solutions violating the analogous condition to the DBG. The following section presents 
several numerical tests of the artihcial forces to speed up convergence to the desired 
distribution. Before closing the present paper, we show the mathematical assurance 
of the faster convergence by analyzing the eigenvalues of the operators in the Fokker- 
Planck equation. In the last section, we summarize our study and provided direction 
for future plans. 

2. Sampling method: Monte Carlo simulation 

Let us briefly review the method to achieve faster convergence to the desired distribution 
in an MGMG simulation from its fundamental ingredients 0. 

2.1. Master equation 

Markov-chain Monte Garlo simulations imitate the discrete-time evolution of the master 
equation given by 

^’(y.i + i) = i:i’(y|x)P(x,(), (1) 

X 

where t is a discrete time and P(x, t) is the instantaneous distribution for the degrees of 
freedom x. The transition matrix P(y|x) dehnes the rule of the stochastic dynamics in 
the master equation. At least, the conservation of the probability (GP) demands that 
the transition matrix should satisfy 

EP(yW = i. (2) 

y 

In addition, the nondiagonal elements of the transition matrix must be non-negative. 

2.2. Balanced condition and its solution 

In MGMG simulations, we numerically perform an iterative update of the instantaneous 
distribution following the master equation. The distribution function in the steady state 

p(ss)(x) 

is dehned by the balance condition 

p‘“'(y) = EP(y|x)P‘“'W. (3) 

X 

By use of the transition matrix satisfying the BG, we obtain the desired distribution 
function, which is set to be P^®®)(x), after sufficient iteration steps. For simplicity, we 
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restrict ourselves to the case aimed at generating the Gibbs-Boltzmann distribution 

p(ss)(x) 

= exTp{—(3E{x) + f3F), where F is the free energy and ft is the inverse 
temperature. 

It is not simple to obtain the nontrivial solution under the BC. The most standard 
way to hnd the solution is to use the detailed balance condition 

F(x|y)p('^^)(y) = P(y|x)p('^‘^)(x) (4) 

for y 7^ X. One can conhrm that the DBG indeed satishes the BG by use of GP. However, 
this is not a unique way to hnd the nontrivial solution under the BG. Recent studies on 
DBG violation have led to new way to construct the transition matrix, as introduced 
below. 

2.3. Violation of the detailed balance condition 

The construction of the transition matrix can be roughly categorized as 1. optimization 
following some rule of the transition matrix neglecting the DBG and construction of a 
duplicate system. 

2.3.1. Suwa-Todo method The former type of method was recently proposed by Suwa 
and Todo [9j. They considered the minimization of the rejection rate, which corresponds 
to the diagonal elements of the transition matrix. In numerical implementation of the 
MGMG method, we generate a candidate as the next conhguration at each update. 
Then we accept the candidate following the probability determined by the nondiagonal 
elements. The average of the rejection rate is thus given by the summation over the 
diagonal elements of the transition matrix. The diagonal elements are given by the 
summation of the nondiagonal elements from GP. To eliminate the diagonal elements, 
the Suwa-Todo method thus allocates the nondiagonal elements while neglecting the 
DBG but satisfying the BG and GP. Let us rewrite the transition matrix into the 
probabilistic flow as 


l/(y|x) = P(y|x)P<“>(x). 

(5) 

Then the BG and GP are given as 


yi/(y|x) = pW(y), 

X 

(6) 

j;:V'(y|x) =P<“>(x). 

V 

(7) 


respectively. The symmetry of the probabilistic flow R(x|y) = R(y|x) reads the DBG. 
The Suwa-Todo method hnds the rejection-free solution violating the symmetry of the 
probabilistic flow. We do not follow their method in the present paper but we remark on 
its efficacy from theoretical aspects. To eliminate the diagonal elements of the transition 
matrix, we may change the nondiagonal elements of the transition matrix. As also 
shown below in a different formulation, the asymmetry in the transition matrix can 
accelerate relaxation to the steady state, as shown by the present authors in Reference 
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fn\ . We thus conclude that part of the efficacy of the Suwa-Todo method comes from 
the asymmetry in the transition matrix, which violates the DBC. However, Suwa and 
Todo simultaneously reduce the rejection rate by decreasing the diagonal part. In the 
present paper, we restrict ourselves to the case in which the diagonal elements are fixed. 
In this sense, the efficacy of the Suwa-Todo method is not yet completely understood. 

2.3.2. Skewed detailed balance condition The latter type is the skewed detailed balance 
condition (SDBC) [10]. We prepare a duplicate system by introducing a replica of the 
original system. For each system, we construct transition matrices Pi(x|y) and P 2 (x|y). 
Then we define the following master equation: 

Pi(y, t + l) = Y, ^i(y|x)Pi(x, t) + Qi{y)Pi{y) - Qi(y)Pi(y), (8) 

X 

where Pj(x, t) is an instantaneous distribution for each system. The swap transition 
from i to i (the counter system of i) is given by Qi{y) and Qi{y) for the inverse manner. 
The balanced condition for the duplicate system holds as 

P‘-~Ky) = E«(y|x)-P“‘p) + Qi(y)P'-\y) - (9) 

X 

We here impose that the steady state shares the same distribution as p(®®)(x). One of 
the generic solutions satisfying the BC is known as the SDBC. The transition matrix 
under the SDBC is given by 

P,(y|x)pl“'(x) = fl(x|y)pl“'(y) (10) 

for y 7 ^ X and the transition between each system should simultaneously satisfy 

^i(y|y) - A(y|y) + Qi(y) - Q2{y) = o. (ii) 

In other words, the BC for each system in the ordinary form as in Equation (I3|) is 
violated. Instead, we demand the BC for the duplicate system, in which the transition 
between each system is taken into account. 


2 . 4 . Ichiki-Ohzeki asymmetric transition matrix 


The violation of the DBC can be recast as the asymmetry of the transition matrix. We 
rescale the transition matrix and the instantaneous distribution function as 

P(x,t) 


P(x, t) = 


W(y|x) = 


Y^P(®®)(x) 


=P(y|x)/p(««^. 


Then the master equation can be rewritten as 
-R(y,«+ 1 ) = y W' (y|x)fl(x, t). 


( 12 ) 


( 13 ) 
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Two of the conditions, namely, the BC and CP, that the master equation must satisfy 


are simply written as 


^W(y|x)^P("")(x) = -^/P^'^Ky), 

X 

(14) 

^^p(ss)(y)W(y|x) =-\/ph«)(x). 

(15) 


y 


respectively. In this expression, the DBG is interpreted as symmetry of the rescaled 
transition matrix as iy(y|x) = iy(x|y). The present authors have proven that the 
asymmetry of the rescaled transition matrix, violation of the DBG, leads to faster 
acceleration of relaxation to the steady state then that driven by using the symmetric 
part of the rescaled transition matrix [12]. 

We emphasize that the essential aspect of speeding up the relaxation is the 
asymmetry of the rescaled transition matrix. To introduce the asymmetry, for example 
in the SDBC, the transition to the duplicate system was considered [lO]. In contrast, 
in the Suwa-Todo method [9], reducing the rejection rate and the nondiagonal elements 
were changed from those under the DBG as a result. 


3. Sampling method: Langevin equation 


3.1. Langevin equation and its corresponding Fokker-Planck equation 

To generate the stochastic dynamics, instead of the master equation, we may use the 
overdamped Langevin dynamics dehned as 

dx = A(x)df + \/^(iW, (16) 

where x is the vector standing for the iV- dimensional degrees of freedom and dx is 
the inhnitesimal change of x during dt. In addition, A(x) is a time-independent force 
vector, T is the temperature, and dW is the Wiener process for the A-dimensional 
degrees of freedom. The force is not necessarily given by —grad£'(x) as in equilibrium. 
The purpose is hereafter to hnd an artihcial force to realize faster convergence to the 
desired distribution. 

For achieving this goal, let us write the corresponding Fokker-Planck equation to 
the above Langevin equation: 
aP(x, f) 


dt 


= —divJ(x, f). 


where J (x, t) is the probabilistic flow dehned as 
J(x, t) = (A(x) — Tgrad) P(x, t). 


(17) 


(18) 


3.2. Divergence-free condition 

When the system is in any steady state, the divergence of the probabilistic how vanishes. 
Therefore, we demand the following equality in the steady state: 

0 =-divjl“>(x), 


(19) 
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where 

= (A(x) - Tgrad) p(®")(x). (20) 

We refer to the above equality as the divergence-free condition. Our problem is to hnd 
a solution J^®®)(x) satisfying the divergence-free condition. The above equality, which 
holds in any steady state, is analogous to the BC in the context of the master equation. 
Hereafter the important thing is to hnd a nontrivial solution under the divergence-free 
condition. The solution is violating the DBG but under the BC in the same spirit as in 
the previous studies for acceleration of the MCMC method PEQ]. 

3.2.1. Equilibrium solution We can immediately hnd a trivial solution of the 
divergence-free condition. When the probabilistic how in a steady state vanishes, the 
divergence-free condition is satished: 

J("^)(x) = 0, (21) 

where the superscript has been changed from “ss” to “eq”. This steady state is in 
particular called the equilibrium state because no current exists. Then the force becomes 

A(x) = —gradP(x). (22) 

This is the trivial solution in terms of our purpose. Indeed, the resulting force does not 
violate the DBC. This fact can be conhrmed through a formulation of the transition 
probability. To violate the DBC, we add an additional term to the trivial force. Below 
we show several nontrivial solutions, which can indeed violate the DBC. 


3.2.2. Vector-field solutions Let us move on from the case without current. We demand 
the probabilistic flow in a steady state to take the following form: 

J(^")(x) = 7 B(x)p(^^)(x), (23) 


where 7 is a control parameter on the degree of violation of the DBC and will be detailed 
below. Owing to the divergence-free condition, the vector held B(x) must satisfy 


N 


0 = 7II 


'a[B] 


k=i 


T 0[xjfc^ 


>(ss)/ 


(24) 


where we have used the fact that the distribution function is set to be cx 

exp(—P(x)/T) and the subscript for the bracket denotes the element of the vector. A 
trivial solution is given by 


[B(x)]fc = exp(P(x)/T - P/T). (25) 

Then the probabilistic how becomes a constant current as 

J(««)(x) = 7 I. (26) 


Here 1 denotes a vector with all elements unity. We obtain the resulting force as 
A(x) = —gradP(x) -|- 7 exp(P(x)/T — P/T)l. 


( 27 ) 
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For practical use of this force, we may rescale 'yexp{—F/T) as 7 . However, we have 
to be careful about the instability in using the exponential term, in particular in the 
low-temperature region. 

Let us hnd an alternative to the above obtained force to remove the instability of 
the exponential term. We set the vector held as 
M ( dE dE \ 

where [xj^v+i = [x]i and [x]o = [xj^v- Then the divergence-free condition holds. The 
resulting force is given as 



M dE f dE dE \ 

^ fc = -^rT + ^ bn- 

9[xjfc d[K\k+lJ 

This force drives the system while satisfying the divergence-free condition but neglecting 
the DBG. In this sense, the resulting force is analogous to the case of using the Suwa- 
Todo method in the context of the MCMC method [9]. This type of solution is not 
unique. We may permute the index j = 1, 2,..., as j = T’(l), ^*(2),..., P{N), where 
F(-) denotes a permutation of N indices. To reach the best performance, we may pursue 
an optimal permutation to realize the fastest convergence. This is beyond our scope in 
the present study. Let us further hnd any other nontrivial solutions. 



3.3. Divergence-free condition for a duplicate system 

As in the case of the SDBC 110], let us duplicate the system. Then we may neglect the 
DBG for each system but hold the BG for the composite system. The analogy to the 
SDBG leads to the duplicate Fokker-Planck equation 
aP(Xi,X2,f) ^ 

-diViJi(xi, X 2 , f), (30) 

^ i=l 

where Xj is the iV-dimensional vector for each system denoted by i and Jj(xi,X 2 ,t) is 
the probabilistic how dehned as 

Jj(xi, X 2 , t) = (Ai(xi, X 2 ) - TgradJ P(xi, X 2 , t). (31) 

The gradient indexed by i is taken for each system. Again we impose the divergence-free 
condition. The divergence-free condition for the duplicate system can be recast as 

0 = -X!divijf"^(xi,x2), (32) 

i=l 

where 

jf"^(xi,X2) = (Ai(xi,X2) -TgradjP^"®^(xi)p("")(x2). (33) 

We impose the steady state of the duplicate system as p(®®)(xi, X 2 ) = P^®®^(xi)p(®®)(x 2 ). 
We can hnd a nontrivial solution by considering the mixture of the forces on the duplicate 
system as 


Ai(xi,X 2) = - gradiP(xi)-F 7grad2P(x2), 
A2(xi,X 2) = - grad2P(x2) - 7gradiP(xi). 


(34) 

(35) 
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Then the probabilistic flow in the steady state is given as 


X 2 ) = 7grad2B(x2)P<“>(xi)P<‘‘>(x2), 

jf>(x„X2) = -7giad,£:(xi)pW(x,)P'“>(x2). 


(36) 

(37) 


One can confirm that the divergence of the above probabilistic flows can vanish 


immediately because the summation becomes zero. 

One may recall the exchange Monte Carlo simulation, in which a multiple system 
with different temperatures is driven simultaneously, while each realization is exchanged 
several times. The technique exhibits extremely outstanding performance with fast 
convergence to the equilibrium state. Naturally, we may expect the existence of an 
analogous Langevin dynamics to the exchange Monte Carlo simulation. We here mention 
this fascinating problem as a subject of a future study, which will be published elsewhere. 

3.4- Path probability and detailed balance condition 

As shown in the previous subsection, we And nontrivial solutions satisfying the analogous 
condition to the BC flT^ and fl32|) . We here confirm the violation of the DBC in 
these solutions. Let us calculate the transition probability during an infinitesimal time 
interval, [t, t + dt]. We obtain 



(38) 


where x[t] is the location at time t. For the duplicated system, a similar computation 
yields the transition probability. We here omit the arguments of the quantities on 
the right-hand side to simplify the expressions. We use the midpoint prescription, 
X = {'x[t + dt] -|-x[t])/2; we take the partial derivative with respect to the location. The 
ratio of the transition probability between the forward and backward processes confirms 
violation of the DBC owing to the existence of the probabilistic flow as follows: 



(39) 


The first term in the exponential function corresponds to the negative of the time- 


derivative housekeeping heat Qhk, while the second term can be regarded as the negative 
of the time-derivative excess heat Qex, i-e., (HUH] 



The existence of the housekeeping heat signals violation of the DBC and, in this sense, 
7 , which is the coefficient in the vector held B, controls the degree of violation of the 
DBC. Also in the case of the duplicate system, we can compute the ratio of the transition 
probabilities for forward and backward processes as 
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Then the time derivatives of the housekeeping heat and excess heat are given as 

Qhk = 7 (x2'^gradiE(x2) - x/grad2^(xi)) , (43) 

Qex = xJgradiE(xj). (44) 

i=l,2 

In the duplicate system, we also conhrm the violation of the DBG from the existence 
of the housekeeping heat. Let us below evaluate the performance of using the artihcial 
force violating the DBG by the above nontrivial solutions. 

4. Numerical tests 

We implement the above artihcial force to introduce the probabilistic how violating the 
DBG. To conhrm that the artihcial force accelerates the relaxation to the steady state, 
we demonstrate simple numerical tests below. We here restrict ourselves to the case 
of the additional force to the duplicate system. Let us consider the system with the 
following double-valley potential energy: 

U{x) = (45) 

The minima of the potential are located at x = —1 and 1. The equilibrium point is x = 0; 
that is, the average over independent runs converges to {x[t]) = W/-^sam in 

t —)• CX3, where x[t] denotes the realization at the time t. We simulate the independent 
^sam = 1000 runs. The initial condition is set to be at x*^^^[0] = 1. To achieve the 
equilibrium point, the particles must climb up the local maximum and reach the other 
minimum of the potential energy. Therefore, it takes a relatively long time to reach the 
equilibrium point. Then we introduce the additional force to accelerate the stochastic 
dynamics and duplicate the system. The time evolution of the Langevin equation is 
evaluated by using the ordinary method known as the Heun scheme [I 6 ]. We set the 
inhnitesimal time as dt = 0.0001. We take an average over independent Wam = 1000 
runs, while taking the time average during At = 0.1. 

Figure [1] shows a comparison between the cases without and with additional force. 
We describe the orbit on the plane of the locations of the duplicate system. The averaged 
locations of the duplicate system are expressed as (xi) and (X 2 ). The case without 
additional force traces the line between ( 0 , 0 ) and ( 1 , 1 ) because any difference does not 
exist in the duplicate system. In contrast, the case with additional force describes a 
different path from the line. The particles seem to wrap around the swelled area of 
the potential energy in the plane of (xi,X 2 ). For weak 7 , the relaxation to the steady 
state is not necessarily faster. However, an increase of the value of 7 leads to faster 
convergence to the steady state than that for the case without additional force. This is 
one of the benehts in violating the DBG, as shown in previous studies in terms of the 
Suwa-Todo method and the SDBG. In addition, we observe the integrated correlation 
time in the steady state. The integrated correlation time is dehned as the collection of 
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the correlation times: 

_ o _ (OiOi+t) — (Oi) (Oi+t) 

n - - ( 02 ) _ ( 0)2 

In the steady state, the dependence of the correlation time on i is lost. In nnmerical 
computation, we cannot sum the correlation time up to an inhnite number. Thus the 
dependence on i remains. We then take an average over i to avoid this dependence. In 
this computation, we use the location as the physical observable. The results are shown 
in Figure [21 We conhrm the decrease of the integrated correlation time. In addition, the 
integrated correlation time reaches a lower limit against any further increase of the value 
of 7 . The decrease of the correlation time gives a reduction of the necessary number to 
efficiently compute the averaged value from the sampling and thus leads to an increase 
in the efficacy of the sampling. This benehcial point is the same as the one found in the 
Suwa-Todo method and the SDBC. 

In addition, the present authors demonstrated a further application of our method 
to a classical spin model, the XY model, which has a critical slowing down in the 
low-temperature region IE]. The proposed method exhibited outstanding performance, 
removing the critical slowing down of the XY model, which often hampers efficient 
computation in the steady state. 

The remaining problem is the mathematical and theoretical understanding of the 
faster convergence to the steady state and the reduction of the correlation time owing to 
violation of the DBG. Fortunately, the former property can be explained by considering 
the eigenvalues of the operator form of the Fokker-Planck equation and transition-rate 
matrix in the master equation. In the next section, we demonstrate the property of the 
operator form of the Fokker-Planck equation in cases without and with additional forces 
and show the shift of the eigenvalues such that it promotes accelerating the relaxation 
to the steady state. 



5. Mathematical assurance: analysis of eigenvalues 

To clarify the origin of the acceleration to the steady state by the additional force, we 
rewrite the Fokker-Planck equation in operator form [TE]. This expression can help us 
to hnd a mathematical understanding of the acceleration of the relaxation to the steady 
state in our technique. 


5.1. Operator form of the Fokker-Planck equation 


First we rescale the instantaneous distribution as 

p(x.«) = Y(5^, 

y/P(ss)(x) 

Then the Fokker-Planck equation in the equilibrium system can be rewritten as 


dP{x,t) 

di 


N 




+ 


1 d'^E 


2 9[> 


-T 


32 

a[x]; 


P(x,«). 


(47) 


(48) 
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Figure 1. Comparison between the cases with and without additional force. From top 
to bottom, the values of 7 take 1, 2, 5, and 10. In the left panels, the red orbit denotes 
the case without additional force. The blue ones represent the case with additional 
force on the plane of xi and X 2 ^ which is the location of a particle in the duplicate 
system. The right panels show the average of the location over independent runs 
over Nsam = 1000. Red plots are for the case without additional force and the blue 
and purple ones stand for each average of the location of the duplicate system with 
additional force. 
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Figure 2. Integrated correlation time. The horizontal axis denotes the value of 7 and 
the vertical one describes the average of the integrated correlation time. The statistical 
error falls within the size of the data points. 


Let us here define the following operators as 

1 


[ajfc = 




2VTd[^]k 

1 dE 


+Vt ® 


Vf 


5[x]fc’ 
d 


(49) 

(50) 


2\/T9[x]fc ''"a[x]fc’ 

where denotes the adjoint of a. These operators satisfy the following commutation 
relation: 


d 




d‘^E 


9[x]fca[x]. 


- k 


and trivial ones as [[a]^ , [a]J = 0 and 
simply expressed in terms of these operators as 

= S'(a^, a)P(x, t), 


J i 


where 


dt 

S'(a^, a) = — a’^'^a. 


(51) 

= 0. The Fokker-Planck equation is 

(52) 


Here we use the transpose of the vectors denoted by (•)"'". The operator on the right- 
hand side in the Fokker-Planck equation is Hermitian. The square root of the Gibbs- 
Boltzmann distribution can be an eigenfunction of this operator with zero eigenvalue. 
By using the same operator expression, we can rewrite the Fokker-Planck equation for 
the nontrivial solution satisfying the divergence-free condition ([29]). The additional force 
modifies the Fokker-Planck equation into 


aP(x, t) 

m 


^(a^, a)P(x, t) 


N 


k=l 



dE 

9[x]fc 


9[x]fc 



P(x,f). (53) 
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From this observation, we find that the effect of the additional force in the Fokker- 
Planck equation can be given by 7 al"'"B(x)P(x, f)/\/T. The nontrivial solution of the 
vector field (l28ll can be expressed by the operators as 


[B(x)]fc = \/T ([a'f]fc-i + [a]fe-i - [a'^'j^+i - [a]fc+i) 


Thus we reach 


= ,v(at,a)F(x,«), 


(54) 

(65) 


where -P(x, f) = P(x, t)/-dP(®®)(x) and 


N 

hF(a^a) = ^(a^a) -7^[at]fc([a];-_i - [a]^+J. 

k=l 

We find the anti-Hermitian part on the right-hand side of the Fokker-Planck equation 
fl5^ . The anti-Hermitian part does not change the steady state because the 
eigenfunction of the original Fokker-Planck equation is that for the annihilation operator 
a. Thus we expect that the existence of the anti-Hermitian part of the operators on the 
right-hand side is closely related to the acceleration of the relaxation. In addition, the 
trivial solution with the exponentially divergent quantity also yields the anti-Hermitian 
part in operator form. 

Furthermore, the above observation is supported from the Fokker-Planck equation 
for the divergence-free duplicate system fl30|l . The Fokker-Planck equation for the 
duplicate system with nontrivial additional force can be rewritten as 


aP(xi,X2,t) 

dt 


= II*^(al,ai)P(xi,X2,f) 

i=l 


N 


-7E 


(- 


dE d 


dE d 


P(xi,X2,f), (56) 


V^[x2]fc 9[xi]fc 9[xi]fc 9[x2]fe y 

where P(xi,X 2 ,f) = P(xi, X 2 , Here we define the operators a* and aj 

for each system. The commutation relation can be replaced with 

d‘^E 


d' 




V 


d[Ki]kd[Ki]i 


(57) 


and trivial ones as 




= 0 and 


d' 


J k 


d' 


J I 


= 0. We here use the Kronecker 


delta 6 ij, which takes unity when i = j and otherwise vanishes. Thus we reach 


aP(xi,X2,t) 

dt 


= H^({aJ},{ai})^(xi,X2,f), 


where 


{a*}) = E a^) - 7 a^ ai - a{ a 2 . 


1=1 


(58) 


(59) 


The artificial force satisfying the divergence-free condition involves the anti-Hermitian 
part also in the case for the duplicate system. From this common property, the anti- 
Hermitian part of the operators in the Fokker-Planck equation is expected to accelerate 
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the relaxation. In addition, when we violate the DBG in the master equation, namely, 
lh(x|y) 7 ^ hh(y|x) in the expression in Equation ([13]), we hnd that the transition 
matrix takes the asymmetric form. The key point comes from the asymmetry of the 
transition matrix in the master equation or operators in the Fokker-Planck equation. 
The relaxation time is characterized by the gap between the hrst and second largest 
eigenvalues of the transition matrix and operators. Below let us analyze the eigenvalues 
for asymmetric operators by identifying both of them for simplicity. Then we prove that 
the gap of the cases with the additional force is always opened in comparison to that 
without it. 


5.2. Shift of the second eigenvalues 


Let us consider the eigenvalues for the asymmetric matrix and operators. We restrict 
ourselves to the case of the nontrivial solution of the vector held in the single system. 
The following analysis can be straightforwardly generalized to the duplicate system case. 
The trivial eigenfunction for the operator W (al, a) is y^P(®®)(x) with zero eigenvalue that 
is the largest. Thus relaxation to steady state P^®®)(x) is assured [18]. The relaxation 
speed is dominated by the gap between the hrst and the second largest eigenvalues. As 
seen above, introduction of the artihcial force induces the anti-Hermitian part of the 
operator iy(al,a), while the Hermitian part is hxed. We here utilize a mathematical 
tool from linear algebra known as the Ostrowski-Taussky inequality |T9| 

|dety4| > det——, (60) 

where A is an arbitrary operator with {A + Al)/2 being positive dehnite. We can show 
that the anti-Hermitian part opens the gap larger by applying the Ostrowski-Taussky 
inequality to the characteristic polynomial det — kT(al,a)^, where / is an identity 

operator. Since the largest eigenvalues for both W (al, a) and (w (al, a) -|- a)^ /2 

are = 0, the Ostrowski-Taussky inequality yields 


det (^XI — W{aL\ a)) 


> det (^XI — S'(a^, a)^ 


(61) 


for arbitrary real A > 0. We show a schematic picture of the behavior of two 
characteristic functions in Figure [3l The root of both of the characteristic functions is 
located at A = 0. In addition. Inequality fl^ ensures that the slope of the characteristic 
function of the asymmetric operator at A = 0 is steeper than that of the corresponding 
symmetric part. To prove the fact that the second largest eigenvalue, the nearest 
root to the origin, is located at the further point by induction of the anti-Hermitian 
operator, let us consider det (^XI — IF(al,a)^ /A, which has the same roots as those of 
the characteristic function except for A = 0. We can prove the lemma that the Hermitian 
and anti-Hermitian operators S\ and G exists satisfying 


det (AJ - W{af a)) 
A 


det (^Sx{a\ a) — G(a'f, a)^ 


(62) 







Mathematical understanding of violation of detailed balance condition 


16 



Figure 3. Real parts of det (A/ — W(afa)) and det (A/ — S'(a'l', a)) . Inequality (l64ll 
claims det (A/ — W{afa)) < det (AI — S{afa)) for < A < 0. 


for arbitrary real A using the appropriate unitary transformation. Since the unitary 
transformation does not change the symmetries and determinants of the operators, we 
reach 


det (^Sx{3i\ a)^ 


det (A/-^(at,a)) 
A 


(63) 


For A satisfying < A < 0, where denotes the second largest eigenvalue for 

the symmetric operator, S'A(al,a) is positive definite. Resorting the Ostrowski-Taussky 
inequality, we find 


det (AJ- hF(at,a)) 
A 


det (A/-^(at,a)) 
- A 


(64) 


This inequality shows that the second largest eigenvalue of the asymmetric operator is 
further from the origin than that of the corresponding symmetric operator. All of the 
above ingredients shows that 

ReA^"^^”) < (65) 


Therefore, the acceleration of convergence by the introduction of asymmetry in the 
transition dynamics is guaranteed in terms of the shift of eigenvalues. 


6. Conclusion 

We investigated several solutions of nontrivial forces to accelerate convergence to the 
desired distribution in Langevin systems. The solutions are roughly categorized into 
two types. The hrst is the addition of an exponentially divergent force to the gradient 
of the potential, i.e., natural force. However, this type of force involves numerical 
instabilities such as overflow. In contrast, the second solution, which takes advantage 
of rotation of the vector field, is free from such numerical instabilities. The latter 
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solution can be simplified by employing the idea in the duplicated system inspired by 
the SDBC. Our proposed method has been indeed confirmed by numerical tests to 
accelerate convergence. This aspect of violation of the DBG can give mathematical 
assurance in terms of the shift of eigenvalues for asymmetric operators. In numerical 
experiments, we also confirm that our method shortens the correlation time compared 
to that in the ordinary implementation of Langevin dynamics. Previous studies on using 
MCMC method without the DBG also revealed the reduction of the correlation time, 
and we may thus regard this property as a beneficial aspect of the DBG violation. This 
advantage might be explained by some mathematical property on the violation of the 
DBG, but we have not reached any conclusive result in this regard. We hope that a 
future study will clarify the reduction of the correlation time for the case violating the 
DBG. 

In the present study, we utilize a mathematical tool in linear algebra to prove 
the acceleration of convergence to the steady state. However, recent developments in 
nonequilibrium statistical mechanics can explain two of the beneficial points of the 
violation of the DBG. The acceleration of convergence and reduction of the correlation 
time rely on the change of frequency of state transitions. Let us consider the case with 
local minima of potential energy. The existence of local minima can be a bottleneck in 
the convergence toward the equilibrium state. This is because the typical events are to 
be trapped in the local minima. As shown in our present study, violation of the DBG 
can help the particles escape from local minima and accelerates the relaxation to the 
steady state. This implies that the particles pass through the modified path from the 
typical transitions in the equilibrium case. The modification of the path probability 
between different microscopic states by additional forces has been recently discussed in 
the context of biased sampling in nonequilibrium statistical mechanics [20l |2T]. The 
detailed analysis on the effect of the additional force will be reported elsewhere from 
the viewpoint of biased sampling and its optimality by considering a kind of variational 
principle. Such a viewpoint will give a deeper understanding on the faster convergence 
to the steady state resulting from the violation of the DBG and fundamental aspects 
on nonequilibrium statistical mechanics through concrete instances near the equilibrium 
system. 
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